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Abstract 

The interpretation of experimental measurements at the LHC requires accurate theoretical predictions for exclusive 
observables, and in particular the summation of soft and collinear radiation to all orders in perturbation theory. We 
report on recent progress towards the automated calculation of multi-parton LHC cross sections at next-to-leading 
order in QCD, including the summation of next-to-leading logarithmic corrections through the combination with 
parton showers. 
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1. Introduction 

Theoretical calculations within fixed-order perturba¬ 
tion theory allow for accurate predictions of inclusive 
observables like total cross sections. The analysis and 
interpretation of experimental signatures at the LHC, 
however, require theoretical predictions for exclusive fi¬ 
nal states, i.e. predictions for differential distributions or 
cross sections with cuts on kinematic variables. Higher- 
order calculations for such exclusive final states in¬ 
volve in general large corrections from soft or collinear 
parton emission, which need to be summed to all or¬ 
ders. An efficient way to achieve such a summation 
is through a parton shower. Parton showers typically 
include the leading logarithmic contributions from soft 
and collinear gluon or quark emission to all orders in 
perturbation theory. They form a central part of Monte 
Carlo event generators and are thus essential to con¬ 
nect theoretical models with realistic experimental sig¬ 
natures. On the other hand, standard parton shower 
event generators often rely on leading-order expressions 
for the hard scattering processes and can therefore not 
predict inclusive cross sections accurately. 

A central goal of recent and current theoretical work 
in LHC phenomenology is thus the automated calcula¬ 
tion of next-to-leading (NLO) LHC cross sections in¬ 
cluding the summation of large corrections from multi¬ 


ple quark and gluon emission through parton showers. 
Such a calculation should combine the accuracy of NLO 
predictions for inclusive cross sections with the power 
of parton shower Monte Carlo programs to reliably de¬ 
scribe differential distributions and cross sections with 
cuts on kinematic variables. A naive combination of 
NLO calculations with parton showers would, however, 
lead to double counting of higher-order contributions 
that are included in both the NLO cross section and the 
parton shower. These contributions have to be identi¬ 
fied and subtracted from the calculation by means of a 
matching procedure. 


Much work in recent years has been devoted to the 
formulation of matching schemes that allow for a con¬ 
sistent combination of parton showers and NLO calcu¬ 
lations including loop-corrections [1-16]. So far, how¬ 
ever, only few publicly available computer codes exist 
that implement such schemes, including in particular 
Mc@NLO [17], Powheg [18], Sherpa [19] and Mad- 
graph5_aMc@NLO [20]. The problem of double count¬ 
ing contributions that are included in both the NLO cal¬ 
culation and the parton shower has been solved in all 
approaches in a similar way: The parton shower evolu¬ 
tion is generically described by a Sudakov factor of the 
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form 

exp (“X.¥X*S /, " 4<?J ' z) )- 

shown here for gluon emission off a quark with virtu¬ 
ality Q 2 . In the collinear limit q 2 —> 0 the function 
P g/t} (q 2 ,z) reduces to the Altarelli-Parisi splitting func¬ 
tion. The Sudakov factor includes short-distance con¬ 
tributions which are also part of the NLO cross section 
and which can be identified by expanding the exponen¬ 
tial in powers of a s . When combining parton showers 
with NLO calculations these contributions need to be 
subtracted to avoid double counting. 

Early work by many authors has focused on imple¬ 
menting specific processes into the Mc@NLO program 
and within the Powheg framework. While the orig¬ 
inal Mc@NLO code is very successful and now in¬ 
cludes a large number of processes, it is tied strongly to 
the parton shower of Herwig (both Herwig6 and Her- 
wig++). This restriction has been lifted recently within 
Madgraph5_aMc@NLO, which is not only fully au¬ 
tomatic, but also parton shower independent and may 
be used with various versions of Herwig and Pythia. 
We note that Sherpa contains an implementation of 
Me@NLO as well. Powheg, on the other hand, allowed 
for matching with a generic parton shower from the be¬ 
ginning and could thus be interfaced with various Monte 
Carlo programs. In practice, the implementation of spe¬ 
cific processes has been greatly simplified thanks to the 
Powheg-Box tool. In each of the frameworks, multiple 
parton emission is summed at the leading-logarithmic 
(LL) level only (although leading-color subleading log¬ 
arithmic effects might be included in parton showers 
based on NLO subtraction schemes, such as the Catani- 
Seymour subtraction), and the results can thus not com¬ 
pete with the accuracy of dedicated resummation calcu¬ 
lations which are routinely performed at next-to-leading 
logarithmic accuracy (NLL). 

The formulation and improvement of parton showers, 
which are in general based on various assumptions and 
approximations, has also been addressed by various au¬ 
thors recently [21-37], Notably, Refs. [23, 26, 27, 32] 
have proposed a parton shower that includes quantum 
interference, spin correlations and sub-leading color ef¬ 
fects. Interference is treated in standard parton showers 
only approximately by means of angular ordering, while 
spin information is generally ignored completely. Spin 
correlations, on the other hand, are crucial for example 
to explore new physics models in cascade decays at the 
LHC. While Refs. [23, 26, 27, 32] describe the theoret¬ 
ical formulation of parton showers with quantum inter¬ 
ference, the first implementation [35] so far is based on 


the standard spin-averaged and leading-color treatment 
of parton splitting only. 

The set-up of the parton shower also has implica¬ 
tions for the matching with NLO calculations. If the 
splitting functions of the parton shower, Pjj(qr,z), and 
the momentum mapping follow closely the definition 
of the subtraction terms used to regularize the soft and 
collinear divergences in the NLO calculation, then the 
matching is simplified considerably. With that in mind, 
various authors have formulated parton showers that are 
based on commonly used subtraction schemes [25, 24]. 

To improve on existing NLO plus parton shower im¬ 
plementations, quantum interference contributions, spin 
correlations and sub-leading color effects in the parton 
shower should be included systematically. Ultimately, 
one would like to develop a fast and numerically ro¬ 
bust code for the automated calculation of multi-parton 
LHC cross sections at next-to-leading order, including 
the summation of next-to-leading logarithmic correc¬ 
tions through combination with parton showers. In this 
contribution we shall describe recent progress towards 
this goal. 

The article is structured as follows. In Section 2 we 
shall describe the construction of a NLO subtraction 
scheme as derived from a NLL parton shower, and its 
implementation into the program package Helac-NLO. 
Section 3 reports on the application of such a scheme for 
the automated calculation of NLO-QCD corrections to 
the production of four bottom quarks at the LHC. Sec¬ 
tion 4 finally discusses progress towards developing and 
implementing a NLL parton shower and the matching 
with Helac-NLO. We conclude in Section 5. 

2. Subtraction schemes for NLO-QCD calculations 

In the following, a new subtraction scheme based on 
a parton shower introduced by Nagy and Soper [23] is 
described. This new subtraction has been implemented 
in the publicly available Helac-Dipoles software [38] 
and has already been tested in the calculation of the 
NLO QCD corrections to pp —> bbbb + X production 
at the LHC (see the next section for more details). Our 
main motivation, however, was to provide a framework 
for a simple matching between a fixed-order calculation 
and the new parton shower. However, before address¬ 
ing this issue, the problem of the integration of subtrac¬ 
tion terms over the unresolved phase space needed to be 
solved. This is the non-trivial part of any subtraction 
scheme. Contrary to the usual practice at NLO, we did 
not perform involved analytic integrations, but rather 
used a semi-numerical approach. After a suitable pa¬ 
rameterization, numerical integrations, inspired by re- 



M. Czakon, M. Kramer, and M. Worek / Nuclear Physics B Proceedings Supplement 00 (2015) 1-22 


3 


cent NNLO methods [39], have been performed. This 
has allowed us to cover both massless and massive cases 
with comparable effort. Let us emphasize that our semi- 
numerical approach distinguishes our work from earlier 
publications on the subject presented e.g. in Refs. [40- 
42], 

Let us start with the inclusive NLO QCD cross sec¬ 
tion for a generic process involving m + 1 final-state 
QCD partons with momenta p„ + pi, —» p\ + • • • + p m+ i, 
which can be written as follows 


C"NLO 


- I 
♦ / 


d<& m &f({p)m) F m 


d<L Wt & R ({p) m+ l) F m +l 


f 

f dx f, 

JO Jm 


({/?}„,) F m 

d<!> m (x) Jl c (x, {p\ m ) F m 


where 

3 ? = |M Bon f, 


R — I A/fReali2 


JC = | At 


y\ v = 2 % [M Born (M'- Lo °P)*| , 

and Al Bom , Al 1 ' Loop , Al Real describe the Born, one-loop 
and real-emission matrix elements, respectively. The in¬ 
tegration measure for the m— and the m +1 -parton phase 
space is denoted by and d<l > m+t , whereas F m and 
F m+ 1 are the jet functions. 

For m well separated hard jets, the Born contribution 
is finite, whereas the virtual and the real-emission terms 
are individually divergent due to the presence of soft and 
collinear singularities. For an infrared-safe definition of 
partonic jets, all soft and collinear divergencies that af¬ 
fect the virtual and real corrections should cancel for the 
inclusive cross section, except for the singularities aris¬ 
ing from the emission of nearly-collinear partons off the 
initial state, which are absorbed into a re-definition of 
the parton distribution functions (PDFs). This absorp¬ 
tion is achieved by introducing suitable collinear coun¬ 
terterms, jf. However, the individual pieces, J[ v and 
JR r , still suffer from soft and collinear divergencies and 
cannot be integrated numerically in four dimensions. To 
solve this problem, local counterterms, that are de¬ 
signed to match the singular structure of the integrand 
in the soft and collinear limits, can be introduced: 


tr nlo - 


r 

Jm 


dQ>m Jt B ({p},n) F m 


d<b m+ 1 \^ R ({p}m+\) F m+ 1 - & D ({p} m+ l) F m ] 

Jm+l 


+ f dx f dO,„(x) S(1 - x)[j{ v ({p) m ) + 

Jo Jm L \ J1 

+ Jf(x, \p\ m ) ] F m . 

They are defined on the (m + l)-parton phase space, 
denoted {/;),„+1 and are subtracted from $\ R and added 
back to after integration over the phase space of 
the unresolved parton. This procedure, called subtrac¬ 
tion method, makes the integrals individually conver¬ 
gent and thus well suited for a Monte Carlo integration. 

The construction of the local counterterms is inspired 
by the well known property of the universal factoriza¬ 
tion of QCD amplitudes in the soft and collinear limits. 
The singular structure of an (m + l)-parton squared am¬ 
plitude for two partons p -, and pj that become collinear 
can be expressed as follows: 


(M({p} m +i)\M({p} m +i)) sing ~ 

{Mm ™)I (vt. ■ Vy) \M({p}m ) )) , 

where \Ai({p} m )) is an amplitude for m on-shell exter¬ 
nal partons, Vy is an operator acting on the spin part 
of the amplitude and {/?},„ describes the reduced m- 
parton kinematics in the limit where partons p, and p t 
become collinear. The structure of the real-emission 
contribution can therefore be reduced to the product 
of a finite Born amplitude squared times a divergent, 
collinear splitting kernel O' associated with the split¬ 
ting Pi -> Pi + Pj : 

W R ({, pUt) * ^ B ({pti ) )®C {ii \p i -p h pj), 

where ® denotes spin correlations. When a parton p t 
becomes soft, on the other hand, factorization can be 
written in the following form: 


(^ / ]({p}/H+l)|2Vl({p} M; +l)) s i n g ~ 

Y5 M ^)\ (T/ ■ TQ . 

k*j 

where T,, are operators acting on the color part of the 
amplitude and {/?}„ is the soft limit of the kinematical 
configuration {p) m +\. In the limit where pj —> 0, one 
can write the real-emission contribution as 

fl R ({p}m+l) ~ 'YjJ{ B {{p} ( m)®S ikj \p i ,p k \Pi,p k ,p j ), 

k*j 

where the factorization is expressed in terms of m soft 
splitting kernels S <k] \ one for each external parton, and 
the symbol ® denotes color correlations. 

Factorization properties described so far dictate gen¬ 
eral rules for constructing local counterterms for a given 
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subtraction method. In the first step, a complete set of 
transformations, which map the original m + 1 partons 
phase space, {p},„+i, into a new one, {/?}„,, that describes 
m on-shell partons, needs to be defined. Subsequently, 
a set of splitting functions D (C) {{p} m , {p} m +i), matching 
the behavior of the soft and collinear kernels in the sin¬ 
gular limits, needs to be worked out. Consequently, the 
local counterterms take the general form: 

N 

& D ({p} m+ 1) = ■ 

(=1 

There is a freedom in defining both mappings 
and splitting functions away from the singular limits, 
each choice leading to a different subtraction scheme. 
The most widespread version is the so-called Catani- 
Seymour scheme (CS) [43, 44], where the local coun¬ 
terterms take the following form: 

3^ cs({f)<»+i) = 

m +1 

2 ® mm k \{p)m + 1 ) 

i.j,k=l 

Here, each mapping \p)\',! k> is labeled by three parton 
indices. For a large number, m , of external partons, the 
number of mappings and matrix elements required by 
the calculation scales cubically: Acs ~ m 3 . For our 
new Nagy-Soper subtraction scheme (NS), on the other 
hand, we can write 

m+ 1 

ij,k= 1 






A 


) 


Therefore, each mapping is characterized by two labels 
{ij } only, which implies that N^s ~ m 2 , i.e. the num¬ 
ber of mappings and subsequent matrix element evalu¬ 
ations is reduced by a factor m compared to the Catani- 
Seymour case. 

Our particular choice of the mapping and a form 
of splitting functions for the Nagy-Soper subtraction 
scheme can be found in Ref. [38]. Let us only mention 
that the latter are based on the original matrix elements 
for q —» qg, g —> qq and g —> gg. Let us also add here 
that the general structure of the real-emission contribu¬ 
tion to cr nlo in this scheme takes the following form: 


ore 


r 

Jm +1 


m +1 

- J] F m ] 

<44=1 

m p 

+ J] d<S> m ft B {{p} m ) ® l <lk) (e, {p } m ) F m 

i,k= 1 

m M ~ 

+ Z Z L dx I d<t> m (x) {p},,,) 

i={a,b ) k= 1 ^ m 

® \K (lk HxAp} m ) + P ( ' k) (x,p 2 F )\F mJ 

where 1(e) and K/P correspond to the integrated sub¬ 
traction terms. More precisely 1(e) encodes the full 
soft/collinear structure of the matrix element in the form 
of single and double poles in e = (d - 4)/2, with d 
the number of space-time dimensions, together with a 
finite part. The K/P operator consists of purely finite 
pieces coming from the initial-state splitting (K) as well 
as from the collinear counterterms (P) and involves an 
additional integration over the momentum fraction x of 
an incoming parton after splitting. 

An important feature of the Nagy-Soper mapping, 
that we would like to emphasize, is a change of all 
spectator momenta at the same time. This fact im¬ 
pacts the way the factorization of the phase space is 
performed. While in the case of initial-state emission 
there are no substantial differences with respect to the 
Catani-Seymour scheme, the factorization is derived in 
a slightly different way when the splitting occurs in the 
final state. The Lorentz invariant phase space for a 
generic final state with m + 1 partons is organized in 
terms of recursive splittings 


d^m+l (.Pi, Pji ) ’ ’ - k n[ — 1. (I) 


d 2 pi d 2 pj 


d 3 k i 


d 2 k ,, 


1-1 


(2n) d 2p° (Inyip] (ln) d lk\ (2n) d 2k° m _ l 
(2n) d 6 d (Q - pi -pj-ki - k m -i) 


r K <"« dK 2 r 
Jk 2 . 2 tt J/i 


dPf. 

I ,k m -r,K) X 

Zn 


d<t> 2 (Fij,K: Q)d0 2 (pi,pf,Pij ), 


where 

p ij=Pi+Pj, Q = K + Pjj, 

and K is the so-called collective spectator momentum 
that is built by all spectator momenta 



d^Pni -1 [ AI {{p} m + 1) F m+ \ 
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Moreover 

^min = ('»ii + ■■■ +»4-,„_l) 2 

k Lx = (VG 2 - m P, - m Pj ) 2 

P iin = (» l P, + m pf 

Pi ax = (VC^-VF) 2 . 

The masses of the on-shell final-state partons are de¬ 
noted by m Pi , m pj , , while Va' 2 is the invariant mass 
of the collective spectator. One can observe that 

d^ni— 1 (&1 > ’ ’ ’ 5 km- 1 5 A") — £/O m _ 1 (&1, * , km-11 A') i 

which follows from the fact that the mapping K —> K is 
a Lorentz transformation, and the phase space is Lorentz 
invariant. In the frame where the total momentum Q is 
at rest, the two-body phase space can be parameterized 
in terms of angular variables, 

1 A(Q 2 ,P 2 K 2 )^ r 

d$>p(P u ,K\ Q) = -—- P-- - , 

u ’ 8 (2n) d ~ 2 (g 2 )¥ J 

where represents the solid angle in d dimensions 

and A is the standard Kallen function, 

A(x,y,z) = x 2 +y 2 + z 2 - 2 xy - 2xz - 2 yz . 


When the total momentum Q is at rest, the in¬ 
tegral fd Q f /_i for the two phase space elements 
d^iiPij, K\ Q) and d<& 2 (pi, K\ Q) is the same. This im¬ 
plies that the Jacobian related to the mapping P t j —> p, 
is given by 


d(t> 2 (Pij,K- Q ) 


A(Q 2 ,Pf r K 2 )^ 
A(Q 2 , m 2 , K 2 ) 


d<A> 2 (pi,K; Q) . 


The phase space for the final-state emission can there¬ 
fore be written in the fully factorized form 


d&m+i(ph Pj. k \, ■ G) — 
d®m(Phk l^’’’ i k / H -], G) x dipf in , 

where 


dPfj 

d^n- 2n 


(A(Q 2 ,P 2 ,K 2 )) 
A(Q 2 ,m 2 , K 2 ) 


d<& 2 (pi,ppP Pi + yQ) 


is the measure of the splitting phase space in d dimen¬ 
sions. The parameters /? and y are uniquely fixed by 
setting 

G = G , K 2 = K 2 , p 2 = m 2 . 


and are given by 


P = 2 


(Pij • G) 2 - p2 j Q 2 


\ (m 2 + 2 P t j ■ Q - P 2 ) 2 - 4 m 2 Q 2 


2 P tJ ■ Q+ P(P 2 j-2 P,j ■ Q-m 2 ) 

2G 2 


Therefore, in the singular limit one would simply have 
P 2 j = m 2 ,p = 1, y — 0 and p t = P,j. 



Figure 1: Parameterization of the angular variables for the final-state 
splitting pi —> pi + pj. Here pk is the spectator parton selected to 
define the azimuthal variable (pj. 

The phase space of the splitting has three degrees of 
freedom in d -4 dimensions. One possible way of 
parameterization is to use Lorentz-covariant scalar pro¬ 
ducts and splitting variables as proposed in Ref. [40]. 
This simple choice results in compact formulae of the 
integrated dipoles for massless partons [41], Applying 
the same strategy to the fully massive case, however, 
the kinematical bounds of the splitting become much 
more complicated and the resulting expressions turn out 
to be very cumbersome. To keep the final expressions 
reasonably compact, we have adopted an alternative pa¬ 
rameterization. For example, in the case of the final- 
state emission, for a set of momenta [p} m the splitting 
Pi —* Pi + Pj and the set of momenta {p} m+ \ had to be 
constructed out of three parameters that we have called 
collinear, soft and azimuthal variables. In the singu¬ 
lar limit, where P 2 . —> m 2 , the collinear and the soft 

'J ' 

variables have to correspond to the relative angle be¬ 
tween the nearly-collinear partons, f),j, and the energy 
of the unresolved parton, Ej. On the other hand, the az¬ 
imuthal variable, (pj, which is the second angular para¬ 
meter uniquely fixes the kinematics of the splitting. Our 
reference frame consists of an orthogonal set of axes 
(x,y,z), where the z-axis is identified with the spatial 
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direction of the vector p,, as shown in Figure 1. The 
azimuthal variable is then the angle, which separates 
the unresolved parton pj from the the x-z plane. Since 
there is complete freedom in selecting the direction of 
the x-axis, we have decided to place it in the plane />,- 
pk, where pk is the momentum of the spectator. In the 
frame, where the total momentum Q is at rest, the defi¬ 
nition of the soft variable is as follows: 

E = __ 

7 " Ff. - mf + 2 Ve^(P? - cos Gj \p,\) ' 

It is convenient to divide this soft variable by its kine¬ 
matically allowed maximum value to get a normalized 
soft variable e: 

e = Ej/Ef “ . 

The integration of the splitting phase space in case of 
the final-state emission runs over the variables 

e e [0,1] , 

c = cos Gj e [-1,1], 

<p = 4>j e [0, 2n\ . 

The soft and collinear limits correspond to e —> 0 and 
c —> I, respectively. 

For the calculation of the integrated subtraction terms 
in the Nagy-Soper scheme, we have adopted the spin- 
averaged version of the splitting functions as described 
in Ref. [26], Our goal is the integration in d = 4 - 2e 
dimensions over the whole phase space of the splitting. 
However, as a consequence of the increased complex¬ 
ity of the mapping, a fully analytic evaluation of the 
integrals turns out to be demanding. Thus, as an alter¬ 
native we have used numerical approaches to integrate 
over the splitting phase space. More precisely, we have 
decided to adopt a semi-numerical approach to consider 
analytic integration when possible, and Monte Carlo in¬ 
tegration otherwise. Since the general dependence of 
the integrands on the azimuthal variable <p was simple 
and all the azimuthal integrals could be classified into 
three groups, we carried out this part of the integration 
analytically. The dependence on the soft and collinear 
variables was not as simple and led to more complicated 
expressions that we treated numerically. Further details 
can be found in the original publication [38], 

2.1. Implementation in Helac-Dipoles 

We have incorporated the new subtraction method 
based on the Nagy-Soper formalism into the Helac- 
Dipoles package, preserving at the same time all the 


optimizations already available in the code. For a de¬ 
tailed description of the package functionalities, we re¬ 
fer to the existing literature [45, 46]. All elements of 
the calculation that do not dependent on a specific sub¬ 
traction scheme, like the Born matrix elements and the 
color correlators, were already provided by the frame¬ 
work of Helac-Dipoles. This fact has dramatically sim¬ 
plified our implementation. 

The construction of the Nagy-Soper subtraction terms 
is dictated by the form of the splitting functions. They 
contain generic spinors and polarization vectors, which 
enables them to treat simultaneously fixed helicities as 
well as random polarization states. We have provided 
random polarization sampling as a further option avail¬ 
able for the Nagy-Soper scheme. This is an alternative 
to the existing random helicity sampling optimization, 
which uses stratified sampling over the different (inco¬ 
herent) helicity assignments of partons [45]. The option 
for the spin sum treatment can be controlled by the user 
in the configuration file dipoles. conf as described in 
the Appendix of Ref. [46]. 

Besides random polarization sampling, which is an 
important speedup in every calculation, we have also 
provided random sampling over color, or color Monte 
Carlo, for the subtracted real radiation part. This func¬ 
tionality has provided an important speedup for matrix 
elements with a large number of colored external states. 
The general ideas from [47, 48] have been adopted, 
which are also an essential ingredient of the Helac- 
1 Loop package [49]. 

Our implementation of the Nagy-Soper subtraction 
scheme (NS) has been tested and compared to the 
Catani-Seymour subtraction scheme (CS) for some spe¬ 
cific processes. More precisely proton-proton collisions 
at the LHC with a center-of-mass energy of 8 TeV 
have been considered and the following partonic sub¬ 
processes gg —> ttbbg, gg —> ttttg, gg —> bbbbg, gg —> 
ttggg have been studied. They give dominant contribu¬ 
tions to the subtracted real emissions at <9(a^) for the 
corresponding processes pp —> ttbb + X, pp —» tttt+X, 
pp —» bbbb + X and pp —> ttjj + X. Moreover, they 
represent a high level of complexity and test almost all 
aspects of the software, as they involve both massive 
and massless states. We have imposed basic selection 
cuts on jets 

Pr(j) > 50 GeV, \y(j)\ < 2.5, A R(jj) > 1 , 

which have been defined through the anti -kj jet algo¬ 
rithm [50] with radius parameter R — 1. The mass of 
the top quark was set to m, — 173.5 GeV and the bottom 
quark was considered to be massless. Results have been 
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presented for the NLO CT10 parton distribution func¬ 
tions [51] with five active flavors and the correspond¬ 
ing two-loop a s . The renormalization and factorization 
scales were set to the scalar sum of the jet transverse 
masses 

h t = m T (ji ), 
i 

where for the top quark 

m T (t) = y]mj + pf(t) 
and for light jets (also tagged bottom-jets) 
m T (j) = PtU) ■ 

A factor of 1 /4 has been included in the scales for all 
but the gg —> bbbbg process where /jr = /jf - Hr has 
been chosen instead. 

In the following, a few examples from the compari¬ 
son of both schemes are given. If not specified other¬ 
wise, full summation over all color configurations has 
been assumed together with random helicity sampling. 
In Table 1, for example, the total number of subtrac- 


Process 

Number of 

Number of 

Number of 


Dipoles (CS) 

Dipoles (NS) 

FD 

gg -> ttbbg 

55 

11 

341 

gg -> ttttg 

30 

6 

682 

gg —> bbbbg 

90 

18 

682 

gg -> ttggg 

75 

15 

1240 


Table 1: Number of Catani-Seymour (CS) and Nagy-Soper (NS) sub¬ 
traction terms for dominant partonic subprocesses contributing to the 
subtracted real emission contributions a! Oioy ) for the pp —> ttbb+X, 
pp —> tilt + X, pp —> bbbb + X and pp —> ttjj + X processes at the 
LHC. The number of Feynman diagrams (FD) corresponding to the 
subprocesses is given as well. 

tion terms that are evaluated in both schemes is shown. 
Also given is the number of Feynman diagrams corre¬ 
sponding to the subprocesses under scrutiny to under¬ 
line their complexity. For each of the 2 —» 5 processes, 
five times less terms are needed in the NS subtraction 
scheme compared to the CS scheme. The difference 
corresponds to the total number of possible spectators, 
which are relevant in the CS case, but not in the NS case. 

Real emission cross sections are presented in Table 2, 
again for the CS dipole subtraction and the new NS 
scheme. All results have been obtained with the same 
Monte Carlo statistics and the resulting relative errors 


Process 

o'rh [Pb] 

tr RH [Pb] 

gg -» ttbbg 

(28.39 ± 0.04) • 10- 3 

(28.59 ± 0.06) • 10~ 3 

gg -» ttttg 

(16.98 ± 0.02) • 10~ 5 

(17.01 ± 0.03) • 10~ 5 

gg -> bbbbg 

(66.24 + 0.16) • 10~ 2 

(66.06 ± 0.22) • 10~ 2 

gg -> ttggg 

(87.96 ± 0.07) • lO^ 1 

(88.16 + 0.08)- IQ’ 1 


Table 2: Real emission cross sections for dominant partonic subpro¬ 
cesses contributing to the subtracted real emissions at 0(a 5 s ) for the 
pp —» ttbb + X, pp —> tttt + X, pp —> bbbb + X and pp —> ttjj + X 
processes at the LHC. Results are shown for two different subtrac¬ 
tion schemes, the Catani-Seymour (CS) dipole subtraction and the 
new Nagy-Soper (NS) scheme, including the numerical error from the 
Monte Carlo integration. 


Process 

t cs [msec] 

t NS [msec] 

/ RE [msec] 

gg -> ttbbg 

24.8 

13.2 

6.5 

gg -» ttttg 

35.7 

18.5 

11.2 

gg -»bbbbg 

26.6 

16.2 

10.1 

gg -> ttggg 

214.8 

108.2 

48.7 


Table 3: The CPU time needed to evaluate the real emission matrix 
element together with all the subtraction terms for one phase space 
point for two subtraction schemes, namely Catani-Seymour, t cs , and 
Nagy-Soper, ? NS . For comparison, we also give the CPU time for the 
pure real emission matrix element calculation, t RE . All numbers have 
been obtained on an Intel 3.40 GHz processor with the Intel Fortran 
compiler using the option -fast. 


Process 

°"rI,col tpb] 

<col [Pb] 

gg -» ttbbg 

(28.35 + 0.14) • 10- 3 

(28.77 + 0.14)- 10~ 3 

gg -» ttttg 

(17.00 + 0.03) • 10~ 5 

(17.01+0.04)- 10~ 5 

gg — > bbbbg 

(65.71 +0.50) • 10~ 2 

(67.00 ± 0.66) • 10- 2 

gg -> ttggg 

(88.04 ± 0.37) • 10- 1 

(87.76 + 0.31)- IQ’ 1 


Table 4: Real emission cross sections for dominant partonic subpro¬ 
cesses contributing to the subtracted real emissions at 0(a 5 s ) for the 
pp —» ttbb + X, pp —» tttt + X, pp —> bbbb + X and pp —> ttjj + X 
processes at the LHC. Results are shown for random color sampling 
for two different subtraction schemes, the Catani-Seymour (CS) dipole 
subtraction and the new Nagy-Soper (NS) scheme, including the nu¬ 
merical error from the Monte Carlo integration. 

are well below 1%. We observe that the difference be¬ 
tween two evaluations of a given cross section is at most 
twice the sum of the corresponding errors. In Table 3, 
the time measured in milliseconds, needed to evalu¬ 
ate the real emission matrix element and the subtrac¬ 
tion terms for one phase space point is shown. The NS 
is scheme is typically twice as fast as the CS scheme, 
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Process 

^RE.POL IP* 3 ! 

88 -» ttbbg 

(28.50 ± 0.06) • 10^ 3 

88 -» ttttg 

(17.01 ± 0.03) • 10^ 5 

gg -»bbbbg 

(66.23 ± 0.20) • 10^ 2 

88 -> ttggg 

(88.16 + 0.07)- 10-' 


Table 5: Real emission cross sections for dominant partonic subpro¬ 
cesses contributing to the subtracted real emissions at 0(a\) for the 
pp —» ttbb+X. pp —> tttt+X. pp —> bbbb+X and pp —> ttjj+X pro¬ 
cesses at the LHC. Results are shown for random polarization sam¬ 
pling for the new Nagy-Soper (NS) subtraction scheme, including the 
numerical error from the Monte Carlo integration. 

but still a factor of about two slower than the evalua¬ 
tion of the real emission matrix element. Overall, both 
schemes, with their different momentum mappings and 
subtraction terms, show a comparable performance and 
give the same results for total real emission cross sec¬ 
tions. 

The performance of Monte Carlo sampling over color 
and polarization has also been studied. In Table 4 we 
present real emission cross sections, which have been 
evaluated with random color sampling, for both subtrac¬ 
tion schemes. We observe agreement with the results 
presented in Table 2, where a summation over all color 
flows has been performed. One should note that in the 
case of the MC summation the absolute errors are 3 - 4 
times higher. In order to obtain the same absolute er¬ 
rors as for results including a summation of color flows, 
9-16 times more events need to be evaluated. How¬ 
ever, the average number of color flows corresponding 
to a random color configuration, which is evaluated per 
phase space point, is dramatically reduced. The over¬ 
all time to obtain the same result is therefore substan¬ 
tially shortened. Our conclusion is thus that random 
color sampling is a powerful approach, especially for 
processes where the number of gluons is higher and ex¬ 
ceeds the number of quarks. 

Finally, in Table 5 real emission cross sections for 
random polarization sampling for our new NS subtrac¬ 
tion scheme are shown. They should be compared to the 
numbers given in Table 2, where we have used random 
helicity sampling. Perfect agreement is found. 

To conclude this section, a complete implementation 
of the Nagy-Soper subtraction scheme both for mas¬ 
sive and massless partons is now available in the Helac- 
Dipoles software. By design, the Nagy-Soper scheme 
has less kinematical mappings and is, therefore, faster. 
On the other hand, we have observed that the absolute 
error of the most costly (in terms of computational time) 


subtracted real emission contribution was slightly worse 
for Nagy-Soper than for Catani-Seymour. In the end, 
we conclude that both schemes are similar in terms of 
efficiency. We did not consider differences below a fac¬ 
tor of two in error or time, which are moreover process 
dependent, a reason to prefer either scheme. There are 
two advantages of our implementation: First, we can 
now perform better tests when calculating fixed order 
NLO QCD corrections by computing real radiation in 
two different schemes. A case study is described in the 
next section. Second, the integrated subtraction terms 
facilitate the matching of the fixed order calculation and 
the Nagy-Soper parton shower with quantum interfer¬ 
ence. This part is described in Section 4. 

3. A case study: NLO-QCD corrections to the pro¬ 
duction of four bottom quarks at the LHC 

The production of four bottom quarks, pp —> bbbb + 
X, is an important background to various Higgs analyses 
and new physics searches at the LHC, including for ex¬ 
ample Higgs-boson pair production in two-Higgs dou¬ 
blet models at large tan p [52], or so-called hidden valley 
scenarios where additional gauge bosons can decay into 
bottom quarks [53]. Accurate theoretical predictions 
for the Standard Model production of multiple bottom 
quarks are thus mandatory to exploit the potential of the 
LHC for new physics searches. Furthermore, the calcu¬ 
lation of the NLO QCD corrections to pp —» bbbb + X 
provides a substantial technical challenge and requires 
the development of efficient techniques, with a high de¬ 
gree of automation. In Ref. [54] we have performed 
an NLO calculation of bbbb production at the LHC 
with the Helac-NLO system [46]. In particular, we 
have presented results based on the Nagy-Soper sub¬ 
traction scheme introduced in Section 2. Two calcula- 
tional schemes have been employed, the so-called four- 
flavor scheme (4FS) with only gluons and light-flavor 
quarks in the proton, where massive bottom quarks are 
produced from gluon splitting at short distances, and 
the five-flavor-scheme (5FS) [55] with massless bottom 
quarks as partons in the proton. At all orders in pertur¬ 
bation theory, the four- and five-flavor schemes are iden¬ 
tical, but the way of ordering the perturbative expansion 
is different, and at any finite order the results do not 
match. Comparing the predictions of the two schemes 
at NLO thus provides a way to assess the theoretical un¬ 
certainty from unknown higher-order corrections, and 
to study the effect of the bottom mass on the inclu¬ 
sive cross section and on differential distributions. First 
NLO results for pp —> bbbb + X in the 5FS have been 
presented in Ref. [56]. In Ref. [54] we have not only 
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provided an independent calculation of this challenging 
process with a different set of methods and tools, but 
also a systematic study of the bottom quark mass ef¬ 
fects by comparing the 5FS and 4FS results. We note 
that NLO results for the production of four top quarks 
in hadron collisions have been discussed in Ref. [57], 

The calculation of the process pp —> bbbb + X at 
NLO QCD comprises the parton processes gg —» bbbb 
and qq —> bbbb at tree-level and including one-loop 
corrections, as well as the tree-level parton processes 
gg —> bbbb + g, qq —> bbbb + g, gq —> bbbb + q and 
gq —» bbbb+q. In the four-flavor scheme q e {u, d, c, s], 
and the bottom quark is treated massive. The bottom 
mass effects are in general suppressed by powers of 
mi,/p, where p is the hard scale of the process, e.g. 
the transverse momentum of a bottom-jet. Potentially 
large logarithmic corrections oc In (nib/p) could arise 
from nearly collinear splitting of initial-state gluons into 
bottom quarks, g —» bb, where the bottom mass acts 
as a regulator of the collinear singularity. This class 
of ln(»!fo//i)-terms can be summed to all orders in per¬ 
turbation theory by introducing bottom parton densities 
in the five-flavor scheme. The 5FS is based on the ap¬ 
proximation that the bottom quarks from the gluon split¬ 
ting are produced at small transverse momentum. How¬ 
ever, in our calculation we have required that all four 
bottom quarks can be experimentally detected, and we 
have thus imposed a lower cut on the bottom transverse 
momentum, prj> > Pr'u- As a result, up to NLO ac¬ 
curacy the potentially large logarithms in the process 
pp —> bbbb+X are replaced by In (mb/p) —> In (p™™/p), 
with nib p n /"b < l- 1 ' an d are thus much less signif¬ 
icant numerically. Therefore, for the process at hand, 
the differences between the 4FS and 5FS calculations 
with massive and massless bottom quarks, respectively, 
should be moderate, but may not be completely negligi¬ 
ble. 

Our calculation has been performed with the au¬ 
tomated Helac-NLO framework [46], which includes 
Helac-Iloop [49] for the evaluation of the numerators of 
the loop integrals and the rational terms, CutTools [58], 
which implements the OPP reduction method [59-62] 
to compute one-loop amplitudes, and OneLoop [63] for 
the evaluation of the scalar integrals. The singulari¬ 
ties for soft and collinear parton emission are treated 
using subtraction schemes as implemented in Helac- 
Dipoles [45], see the discussion in Section 2. The 
phase space integration is performed with the help of 
the Monte Carlo generators Helac-Phegas [64-66] and 
Kaleu [67], including Parni [68] for the importance 
sampling. 


The Helac-Dipoles package has been based on the 
standard Catani-Seymour dipole subtraction formalism 
[43, 44]. We have extended Helac-Dipoles by imple¬ 
menting the new subtraction scheme [40, 41] based on 
the momentum mapping and the splitting functions de¬ 
rived in the context of an improved parton shower for¬ 
mulation by Nagy and Soper [23], as described in Sec¬ 
tion 2. The results presented in Ref. [54] have been the 
first application of the Nagy-Soper subtraction scheme 
for a 2 —» 4 scattering process with massive and mass¬ 
less fermions. 

Below we shall present a number of selected numer¬ 
ical results for the pp —> bbbb + X cross section at the 
LHC at the centre-of-mass energy of = 14 TeV. We 
discuss the impact of the NLO-QCD corrections, and 
study the dependence of the results on the bottom quark 
mass. 

Let us first specify the input parameters and the 
acceptance cuts we impose. The top quark mass, 
which appears in the loop corrections, is set to m, — 
173.5 GeV [69], We combine collinear final-state par- 
tons with pseudo-rapidity |?/| < 5 into jets according 
to the anti-L/ algorithm [50] with separation R - 0.4. 
The bottom-jets have to pass the transverse momentum 
and rapidity cuts prj, > 30 GeV and |v/,[ < 2.5, respec¬ 
tively. The renormalisation and factorisation scales are 
set to the scalar sum of the bottom-jet transverse masses, 
Pr = p f — po — Hr, with 

H t = niRb + ni T b + niR' b + 
and the transverse mass 

m r,b = ym 2 h + p 2 r b . 

For the five-flavor scheme calculation with massless 
bottom quarks the transverse mass equals the transverse 
momentum, m T b — PR,b■ Note that the implementa¬ 
tion of a dynamical scale requires a certain amount of 
care, as the subtraction terms for real radiation have to 
be evaluated with a different kinematic al configuration 
specified by the momentum mapping of the subtrac¬ 
tion scheme. Comparing the results as obtained with 
the Catani-Seymour subtraction and the Nagy-Soper 
scheme, which is based on a different momentum map¬ 
ping, provides an important and highly non-trivial inter¬ 
nal check of the calculation. 

3.1. Massless bottom quarks within the five-flavor 
scheme 

The NLO predictions for the inclusive cross section 
are presented in Table 6 for the NLO MSTW2008 [70] 
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parton distribution function (pdf), with five active fla¬ 
vors and the corresponding two-loop a s . To study the 
impact of the higher-order corrections, we also show 
leading-order results obtained using the MSTW2008 
LO pdf sets and one-loop running for a s . 

Varying the renormalisation and factorisation scales 
simultaneously about the central scale by a factor of 
two, we find a residual scale uncertainty of approxi¬ 
mately 30% at NLO, a reduction by about a factor of 
two compared to LO. The AT-factor, K = cr NL0 /cr L o 
1.37, is sizeable. Note, however, that the /('-factor is an 
unphysical quantity and depends strongly on both the 
default choice of scale and the pdf set [54], 

In Ref. [54] we have also presented predictions for se¬ 
lected differential distributions which are an important 
input for the experimental analyses and the interpreta¬ 
tion of the experimental data. Figure 2 shows LO and 
NLO predictions for the transverse momentum of the 
hardest bottom jet. We also show the theoretical uncer¬ 
tainty through scale variation and the AT-factor as a func¬ 
tion of the transverse momentum. It is evident from Fig¬ 
ure 2 that the NLO corrections significantly reduce the 
theoretical uncertainty of the differential distributions, 
and that the size of the higher-order effects depends on 
the kinematics. For an accurate description of exclusive 
observables and differential distributions it is thus not 
sufficient to rescale the LO prediction with an inclusive 
Af-factor. 


3.2. Massive bottom quarks within the four-flavor 
scheme 

Within the four-flavor scheme, bottom quarks are 
treated massive and are not included in the parton 
distribution functions of the proton. We define the 
bottom quark mass in the on-shell scheme and use 
nib = 4.75 GeV, consistent with the choice made in the 
MSTW2008 four-flavor pdf [71]. 

The central cross section predictions in LO and NLO 
for p = H T using the 4FS MSTW2008 [71] pdf are 
shown in Table 6, in comparison with the 5FS results. 
We observe that the bottom mass effects decrease the 
cross section prediction by 18% at LO and 16% at NLO. 
The residual scale dependence at NLO is approximately 
30%, similar to the 5FS calculation. 

The difference between the massless 5FS and the 
massive 4FS calculations has two origins. First, there 
are genuine bottom mass effects, the size of which de¬ 
pends sensitively on the transverse momentum cut. For 
pi™ = 30 GeV we find a 10% difference between the 



Figure 2: Differential cross section for pp —> bbbb + X at the LHC 
( ffs = 14 TeV) in the 5FS as a function of the transverse momentum 
of the hardest bottom jet. The dash-dotted (blue) curve corresponds to 
the LO and the solid (red) curve to the NLO result. The scale choice 
is pr = pf = po = Hj. The hashed area represents the scale un¬ 
certainty, and the lower panels display the differential K factor. The 
cross sections are evaluated with the MSTW2008 pdf sets. 



PT b [GeV] 

Figure 3: Differential cross section for pp —> bbbb + X at the LHC 
( sfs = 14 TeV) in the 4FS and 5FS as a function of the transverse 
momentum of the hardest bottom jet, normalised to the corresponding 
inclusive cross section. The scale choice is pr = pf = po = Hj, the 
cross sections are evaluated with the 5FS and 4FS MSTW2008 pdf 
sets, respectively. 

5FS and 4FS from non-singular bottom-mass depen¬ 
dent terms. This difference decreases to about 1% for 
pUjl'P = 100 GeV. Second, the two calculations involve 
different pdf sets and different corresponding a s . While 
a 4FS pdf has, in general, a larger gluon flux than a 5FS 
pdf, as there is no g —> bb splitting, the corresponding 
four-flavor a s is smaller than for five active flavors. For 
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pp — > bbbb + X 

o-lo [pb] 

<T NLO [Pb] 

K = cr NLO /cr LO 

5FS 

qq q+58.7 (59%) 

' y -34.9 (35%) 

1 ^y+38.8(28%) 

1 JO. /_30.9 (23%) 

1.37 

4FS 

o a c+49.7(59%) 
-29.6(35%) 

, , o 0+33.3(28%) 
110 -29.0(24%) 

1.40 


Table 6: 5FS and 4FS LO/NLO cross sections for pp —> bbbb + X at the LHC ( sfs = 14 TeV). The renormalisation and factorisation scales have 
been set to the central value po = Hj, and the uncertainty is estimated by varying both scales simultaneously by a factor two about the central 
scale. Results are shown for the 5FS/4FS MSTW2008LO/NLO pdf sets. 


pp —> bbbb + X the difference in a s is prevailing and 
results in a further reduction of the 4FS cross section 
prediction by about 5%. This latter difference should 
be viewed as a scheme dependence rather than a bottom 
mass effect. 

In Figure 3 we present the differential distribution in 
the transverse momentum of the hardest bottom jet, as 
calculated in the 5FS with massless bottom quarks and 
in the 4FS with mi, = 4.75 GeV, normalised to the cor¬ 
responding inclusive cross section. We find that the dif¬ 
ference in the shape of the distributions in the 5FS and 
the 4FS is very small. 

To conclude this section, we have presented selected 
results for the differential cross-sections for pp —» 
bbbb + X at the LHC at the centre-of-mass energy of 
yfs = 14 TeV [54]. We find that the higher-order cor¬ 
rections significantly reduce the scale dependence, with 
a residual theoretical uncertainty of about 30% at NLO. 
The impact of the bottom quark mass is moderate for 
the cross section normalisation and negligible for the 
shape of distributions. The fully differential NLO cross 
section calculation for the process pp —» bbbb + X 
presented in Ref. [54] provides an important input for 
the experimental analyses and the interpretation of new 
physics searches at the LHC. 


4. Parton shower with quantum interference and 
matching 


In this section, we will discuss the Nagy-Soper 
shower in more detail, including its particular imple¬ 
mentation in the C++ library Deductor [35]. Among 
other topics, we will elaborate on the theoretical frame¬ 
work necessary to include quantum interference effects. 
Furthermore, we will point out the inherent ambiguities 
of the approach. Ultimately, we will discuss the match¬ 
ing to fixed order calculations at the next-to-leading or¬ 
der in QCD, and present some results for a non-trivial 
process: the production of a top-anti-top-quark pair in 
association with a jet in hadronic collisions. This sec¬ 
tion is based on [72]. 


4.1. Basic concepts 

We start from a generic 2 —> m process, which is 
defined by two initial state partons a and b and 1, ...,m 
final state particles. Each particle is described by a set 
of quantum numbers to define the flavor /], spin ,v, and 
color c, of the particle and its momentum p,. The ini¬ 
tial state parton kinematics is described by momentum 
fractions, tj a and with respect to the original collid¬ 
ing hadrons, instead of their momenta. Thus, a complete 
parton ensemble can be described by 1 


\p,f, S,c} m = {[tja,-fa, S a ,Ca], [t]b, ~fb, V>,Q>], 

\.P 1 s f\ i 5 1 , C 1 ]> I Pm•. fm * S m , C m ]} . 


The state of the parton shower evolution is described 
by a quantum density matrix p, which gives the ’prob¬ 
ability’ 2 to find a certain parton ensemble {p,f, s, c},„. 
The expectation value of an observable F for any final 
state multiplicity is given by 


o-[f] = 



[d{p,f ),„] 


f a (Pa,p 2 F )fb(Tl h ,p 2 F ) 
4n c (a)n c (b) x flux 


X < M({p , f } m )| F({p, f } m ) | M([p, /)„,)> 

= V —7 f [d{p,f} m ]Tr\p({p,f} m )F({p,f} m )] . 
*—1 ml J 


Here the 1/m! is the symmetry factor for identical par¬ 
ticles in the final state and [d{p,f} m ] is the m particle 
phase space measure. f a /b(n,p 2 F ) are the parton density 
functions evaluated at the momentum fraction // and fac¬ 
torization scale i u 2 . The factor 4 in the denominator is 
the spin averaging factor. n, (i) represents the averaging 
color factor for the initial state partons, where n c (q) = 3 
and n c (g) = 8. The expectation value of F and the trace 
in the second line are meant to be a summation over 


The minus sign for the initial state flavor is just a convention, 
because all partons are considered as outgoing. 

2 Since p is related to color ordered amplitudes, it may become neg¬ 
ative for subleading color configurations. Thus, one cannot naively 
interpret p as a probability. Nevertheless, we use the terminology of 
statistical mechanics. 
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indices in color ® spin space. The quantum density p 
introduced in the second line, is thus given by 

fa(ria,Lt 2 F )Mrib,V 2 F ) 

P({P, J m) ~ 4/Ic ( a)nc(fe) x flux 

x\M({p,f} m ))(M({p,f} m )\ (1) 

|{s. c} m )p({p, /, s', c , s, c} m ) ({s', c '} m I . 

S,C s',c' 

It is useful to define basis vectors \{p,f,s',c',s,c} m ), 
such that 

p({p, f, s',c', s, c} m ) = ({p, f, s',c', s, c},„\p ) . 

One can then write the expectation value of the observ¬ 
able F as 

cr[F] = (F\p) . 

A particularly important observable for the defintion of 
the parton shower is the total cross section measurement 
function. It is defined as 

(1| {p,f,s',c',s,c} m ) = ({s'} m \{s} m ) {{c'} m \{c} m ) . 

The shower evolution equation describes the propa¬ 
gation of the quantum density matrix from some initial 
shower “time”, to, which represents the hard interaction, 
to the final “time” t F in the low energy regime. The fi¬ 
nal shower time t F characterizes the physical scale at 
which parton emissions cannot be described perturba- 
tively anymore. The definition of shower time t is not 
unique and is explained later. The parton shower evolu¬ 
tion will transform a few partons at the matrix element 
level, to a realistic final state with jets typically made of 
many partons. After this evolution, a phenomenologi¬ 
cal hadronization model must be applied. The pertur¬ 
bative evolution itself is described by a unitary operator 
U(t F , to). The observable F, after showering, has the 
expectation value 

cr[F] = (F\p(t F )) = (F\U(t F ,to)\p(t Q )) . 

The unitarity of the evolution operator is a consequence 
of the requirement that it should not change the total 
cross section. Thus (l\U(t F , to)\p(to)) = (l|p(fo))- The 
evolution operator can be obtained from a differential 
equation involving two operators 'Hitt) and 'V(t), corre¬ 
sponding to the concepts of real and virtual corrections 
respectively: 

dU(t,to) 

, = VHi(i) - *vm U(t, to) . (2) 

at 

Here Hi(t) describes the emission of a resolved parti¬ 
cle, i.e. the momenta, flavor, spins and color configura¬ 
tion will change after its application. 'V(f) describes the 


unresolved emission and therefore does not alter mo¬ 
mentum or flavor configurations. Nevertheless it can 
change color configurations, which will affect further 
emissions. For convenience of calculations, the vir¬ 
tual operator can be further decomposed into 'V(t) — 
’VeH) + *Vs(0, where r V F (t) is diagonal in color space, 
while ''Vs (t) is not. Interestingly, the evolution equation 
takes the same form as the time evolution of a statistical 
ensemble in Liouville space 

^ = l -\p(t),H] = Lp(t) , 
ot a 

where the Liouville operator can be identified as L - 
vH,(t)-^m 

Traditional parton showers are constructed using the 
large N c limit. Thus, the state is always color diagonal 
implying 'Vs (t) —> 0. In this case, Eq. (2) yields 

U(t, to) = N(t, to) + I" dr U(t, tYHi(t)N(t, to) , 

Jlo 

with the Sudakov form factor 

N(t, t 0 ) = T exp j- dr 'V(t)| . 

Since 'V(t) is diagonal in the traditional approach, 
N(t, to) is a number and not a matrix in color space. In 
the general case with non-diagonal 'V(t), it is not practi¬ 
cal to exponentiate a matrix in color space. The idea is, 
therefore, to exponentiate only the diagonal color part 
'V F {t), and treat 'V s ( t) iteratively as a perturbation. This 
can be justified by noting that the off diagonal color 
contributions are always suppressed by a relative fac¬ 
tor of 1 /Nl compared to the leading color contributions. 
Therefore, the solution of Eq. (2) with full color evolu¬ 
tion, using the decomposition of 'V(t) — 'Ve(I) + 'Vj(f), 
is given by 

U(t,to) = N(t,to) 

+ f dTU(t,T)[‘H I (T)-'V S (T)]N(T,to). 

J to 

with 

N(t, t Q ) = T exp dr *V £ (t) j . 

4.2. Real and virtual evolution operators 

The real evolution operator, 'Hi, describes the transi¬ 
tion from an m-particle ensemble to an (m + l)-particle 
ensemble. This is achieved by splitting a chosen parton 
into two, which would physically correspond to a decay 
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of a slightly ofF-shell parton. The splitting is constrained 
by flavor conservation, //—>//+ f m + 1 , and momen¬ 
tum conservation, p\ —» pi + p m +\. The description of 
the transition is ambiguous, as only the singular limits 
of amplitudes are uniquely determined in QCD. After 
emitting a particle, it is necessary to correct the mo¬ 
menta in the event in order to ensure momentum con¬ 
servation and preserve the on-shellness of all particles. 
This is done by certain momentum mapping operators, 
which define momenta and flavors of the new ensemble 

[p, fU+i = Ri({p, fU) ■ (3) 

where l e {a,b,\, ...,m). In Deductor, a global mo¬ 
mentum mapping has been chosen. Whenever a par¬ 
ticle is emitted, the momentum of all final state par¬ 
ticles is affected, see the discussion in Section 2. In 
contrast, e.g. the Sherpa [19] parton shower is based on 
Catani-Seymour dipoles [43], which have local momen¬ 
tum mappings. In this case a single parton momentum 
is modified. An explicit description of the original mo¬ 
mentum mapping used in Deductor can be found ei¬ 
ther in Ref. [23] or in Ref. [38], Recently, the initial 
state momentum mapping has been slightly modified. 
A study [31] showed, that the generated px spectrum in 
pp —> Z strongly depends on the momentum mapping 
for initial state parton splittings. For this reason, De¬ 
ductor uses a momentum mapping, which allows for an 
improved resummation of higher-order corrections [36]. 

Besides momentum and flavor mapping operators, 
'Hi(t) contains splitting functions which correspond to 
the factorisation of amplitudes in the soft, collinear and 
soft-collinear limits. In these limits the amplitude can 
be written as 

\M({p,f} m+ i)) = v({p,f} m+ i)\M({p,f} m )) . 

The operator v(\p, f\ m +\) acts in color and spin space. 
The behaviour of amplitudes in singular limits translates 
into a similar behaviour of the density matrix, which can 
be written more explicitly as 

P({pj)m+ 1 ) ~ 

^ T Jifi -» fi + fm+iWj({p,f} m +i)p({p,f} m ) 

l,k 

x V k ({p,f} m+l )T k (f k -> f k + f m+1 ), 

(4) 

where Tj(fi —* fi + f m + 1 ) is an operator in color 
space, while vj({p,f} m + 1 ) is an splitting operator in 
spin space. The general prescription to obtain De- 
ductor’s splitting functions V) has been summarised 


in [38], whereas the complete set of splitting functions 
can be found in [23], 

The approximation in Eq. (4) can be cast into an op¬ 
erator equation, \p m+ \) = £/*S/|p m )- The operator Si 
describes all possible splittings of the emitter parton l. 
'Hiit) is then defined by the splitting operators S/ at a 
fixed shower time 7~/({p, /},„ 

'H I (t) = J]s l S(t-T l ({pJ} m )) , (5) 

i 

where the sum runs over all possible emitters l. The 
shower time T)(\p, /)„,) corresponds to an infrared sen¬ 
sitive scale, discussed in Section 4.4. 

The virtual evolution operator, 'V(f), represents the 
unresolved virtual corrections. Nevertheless its con¬ 
tent is fixed due to the unitarity condition of the parton 
shower. By applying (1| from the left and |p) from the 
right to Eq. (2) one obtains 

{\\<H I {i)-'V(t)\p{t)) = Q, 

which should be valid for any |p(f))- Parton shower 
unitarity corresponds to simplified virtual corrections, 
whose main function is to cancel the divergences of the 
real corrections. In order to obtain an expression for 
'V(t), we write 

(\mt)\{p,f,c\c, s', s\ m ) = (l| < H/0)|{p,/,c , ,c,/,«),„). 

This equation has an ambiguous solution in color space. 
We shall not reproduce here the explicit form for 'V(t) 
used in Deductor. It can be found in [23]. 

4.3. Logarithmic accuracy 

In the previous section, we discussed the real and vir¬ 
tual evolution operators. While they allow for an ex¬ 
act treatment of color, the current implementation of 
Deductor is based on the so-called LC+ approxima¬ 
tion [32], which amounts to only allowing for color 
evolution without non-local modifications of the color 
state. In this case, we identify 'V/Tfj = r V LC+ (t) and 
Ts(0 = A'V(f) = 'V(f) - 'V LC+ 0). 

The color diagonal part of the virtual operator will be 
exponentiated, which yields a parton shower formula¬ 
tion similar to the traditional one. The evolution equa¬ 
tion is given by 

U LC+ (t , t 0 ) = N LC+ (t, t Q )+ 

f dr U LC+ (t, t)'K / lc+ (t)A lc+ (t, t 0 ) , 

Jto 
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where 


N LC+ (t, t 0 ) = exp dT 'y LC+ (r) 

The color off-diagonal splittings are included perturba- 
tively 

U(tJo) = U LC+ (t,t Q )+ 

f dr U(t,T) [A7T/(t) - AT'(r)] U LC+ (T,t 0 ), (6) 

Jto 

where A'H,(t) = 'H,(t) - Pl) x:+ (t) and A<V(t) = 'V(t) - 

ry LC +( f ) 

Let us investigate the logarithmic accuracy of observ¬ 
ables. We suppose that the parton shower is used to cal¬ 
culate an observable O which contains large logarithms 
L of some invariant. For definiteness, we could image 
this invariant to be the transverse momentum of a gauge 
boson. Then (O) has the form 

(O) = Yj c O’ 2 n)<L 2 " + Y c 0h 2n - 1 )a n s L 2n ~ l + ■■■ . 

n n 

We further suppose that a shower with full color, gen¬ 
erated by U{t,to) reproduces all coefficients c(n,2n) 
and c(n,2n - 1) correctly. Then U LC+ (t,to) will re¬ 
produce c(n, 2n) exactly, because the LC+ approxima¬ 
tion is exact with respect to the soft-collinear singu¬ 
larities. One insertion of [ A77 /(t) - A'LYr)] generates 
a contribution ~ a s L because it contains a correction 
for soft wide-angle gluon emission. This term mul¬ 
tiplies contributions of order a" 1 L 2 "* 2 and, therefore, 
corrects the coefficient c(n, 2n - 1). A second insertion 
of [A7T/(r) - AW'fr)] would only affect the coefficients 
c(n, j < 2n - 2). Therefore, one insertion is sufficient 
to obtain NLL accuracy [32], (An illustration of the 
logarithm counting is given in Figure 4.) This implies 
that NLL accuracy is in fact obtained with the evolution 
equation 

t/(f,f 0 ) = t/ LC+ (f,f 0 )+ 

dT U LC+ (t, t) [A7T/(t) - A'V(t)] U lc+ (t, t Q ) . 



4.4. Shower time 

Emissions generated by a parton shower are strongly 
ordered in some kinematic variable in order to correctly 
resum the leading logarithms, ft turns out that the choice 
of the ordering variable, which we call shower time, is 
ambiguous. The essential approximation made in the 



Figure 4: Illustration of the logarithm counting. One step in the ver¬ 
tical direction is given by an insertion of U LC+ and a diagonal step is 
given by an insertion of - A*V(f). One can see, that two inser¬ 

tions of A 44i(t) - A'V(f) only contribute to the coefficient of a"L 2 "~ 2 . 


parton shower description is that, in each step of the evo¬ 
lution, all partons are on-shell. Thus, the parton shower 
time should allow us to neglect the virtuality of a split¬ 
ting parton. Here, we summarize the discussion from 
Ref. [36] for final state radiation. The discussion for ini¬ 
tial state radiation is analogous. 

Consider a mother parton 0 which splits in two par- 
tons, po —> pi+pi- The momentum po is given in light- 
cone 3 variables (+, -,_L), where po denotes the trans¬ 
verse momentum. We denote the virtuality of the parton 
as Vg. The momentum takes the form 


P o 


P, 


Po + m l + v o 
2 P 



The daughter momenta are given in the Sudakov param- 
eterisation as 


Pi 

P2 


zP, 


P\ + m\ - 
2 zP 


->p l 


(1 -z)P 


Pj + m- 2 + vj 
2(1 -Z)P 



From momentum conservation we obtain 


2 ((1 -z)p\ -zpi) 

V n = 


z( 1 -z) 


1 - Z 


1 - z 


In order to be allowed to neglect vf and v, at each step 
of the evolution, we must require 


— «c vj 
z 


1 -z 


«: Va 


3 k ± = {k° ± k 3 )/ V2 and k 2 = 2 k + k~ ^£ 2 
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Inserting the momentum fraction z 

Pi ' Qo 
Po ■ Qo ’ 

, (Po - Pi) ■ Qo P 2 ■ Qo 
Po ■ Qo po ■ Qo 

where the last approximation is valid in the singular 
limit po ~ p\ + p 2 - Here, Qo is the total final state 
momentum. We finally arrive at the conditions 

7 2 2 2 

vf vt, vf vi, 

-«c - and --— «c - . 

2/? i • Qo 2po ■ Qo 2pi ■ Qo 2po ■ Qo 

They should always be fulfilled. Therefore, we define 


A 2 I(Pt ± Pm+lf ~ m,\ 2 

A? =- -Q: 


2 pi ■ Qo 


o 


and enforce the emissions to be ordered in A j. This can 
be achieved by defining the dimensionless shower time 
as 


Ti({p,f} m ) = - log 


(N) 

\Ql) 


Generally, one can choose other ordering variables. E.g. 
Pythia 8 [73] uses the transverse momentum, p T , in the 
parton splittings to order the emissions, Herwig [74, 75] 
uses angualar ordering, and Pythia 6 [76] uses a virtu¬ 
ality ordering. A consequence of the A 2 ordering is an 
enlarged phase space for initial state splittings as com¬ 
pared to p T ordering [36], 


4.5. Ambiguities of the parton shower definition 

As we have already pointed out at various places, the 
construction of the parton shower is not uniquely de¬ 
fined. In the following we list the main ambiguities. 
They should be kept in mind, since future findings might 
require modifications of their solutions. 

Momentum Mappings: The way the parton shower 
distributes momentum among the particles after 
splitting may influence the resummation accuracy. 
In case of Drell-Yan Z-boson production, a study 
[31] showed that a different momentum mapping 
[36] generates a different pj spectrum of the Z- 
Boson. 

Splitting functions: Splitting functions are only re¬ 
quired to reproduce the singular limits of QCD ma¬ 
trix elements, but can have an arbitrary finite re¬ 
mainder. 


Soft partition function: Soft emissions are spread on 
different collinear emitters by means of a partition 
function. The latter is, however, arbitrary. In [27] 
it was argued, that the partition function should not 
depend on the emitter’s energy. A choice, which 
has not been explored, would be to make the parti¬ 
tion function spin dependent. 

Color treatment: In Deductor, the color diagonal part 
of the evolution is exponentiated, whereas the off- 
diagonal part is treated perturbatively. The sepa¬ 
ration depends on the representation of the color 
algebra, if the perturbative insertion of the off- 
diagonal color operator is truncated. In Ref. [77] 
a different approach for a full color treatment in a 
parton shower has been proposed. 

Spin treatment: The spin basis (not discussed in this 
proceedings) is arbitrary. A different choice could 
modify the spin weights in the spin evolution of the 
shower. 

Shower time: The parton shower emissions are 
strongly ordered. As found in Ref. [36], A 2 
ordering provides a wider phase space for initial 
state splittings than ordering in pj. In many cases, 
nevertheless, different ordering variables give 
the same results. A counter example would be, 
e.g. angular ordering vs. pj ordering. Angular 
ordering preserves color coherence of soft gluon 
emissions, whereas pj ordering cannot account 
for this effect [78], 

PDF evolution: Traditional parton showers interface 
LHAPDF [79] to obtain the ratio of PDFs occur¬ 
ring in the backward evolution for initial state ra¬ 
diation. Therefore, those parton showers explicitly 
depend on the PDF kernels they use. Deductor 
tries to minimize this dependence by evolving the 
PDFs according to the shower splitting functions. 
Since the splitting functions are not uniquely de¬ 
fined, the PDF evolution is also not unique. Ad¬ 
ditionally, if massive quarks are assumed in the 
initial state, then the mass dependent terms of the 
splitting kernel depend on the definition of the 
shower time. 

4.6. Matching at next-to-leading order 

Matching NLO calculations with parton showers is 
a widely explored subject and there already exist sev¬ 
eral matching schemes [1-16]. The most popular ones 
are the Powheg method [9, 13] and the Mc@NFO for¬ 
malism [4, 7]. A general comparison between those two 
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major schemes can be found in [16]. We choose to work 
in analogy to Mc@NLO instead of Powheg, because we 
are looking for a general solution, which can be eas¬ 
ily automated. Before we discuss the problems of the 
Mc@NLO formalism and their solutions, we want to 
give a brief overview of the general objectives of parton 
shower matching. This section presents original results 
obtained in [72], 

Independently of the accuracy of the matching we ob¬ 
tain the following benefits: 

Connection to low energy physics: Inclusive distri¬ 
butions are not affected by showering. Neverthe¬ 
less, the evolution of partons down to a scale tp 
allows to include decays of unstable particles, as 
well as the non-perturbative hadronization and 
multiple interactions models. 

Logarithmic accuracy: Infrared sensitive observ¬ 
ables, which are ill defined at fixed order, are 
replaced by finite predictions due to resummation 
of large logarithms generated by collinear, soft 
and soft-collinear splittings. 

Matching at the next-to-leading order gives us further¬ 
more: 

Cross section normalization at NLO: When consid¬ 
ering an inclusive observable F we want to keep 
the fixed order normalization of the cross section. 
Thus, the parton shower must not modify the total 
cross section 

(C|t/(f F ,f 0 )lp(fo)) = cr NLO [F] . 

High-px emission according to matrix elements: 

The parton shower is valid in the soft and collinear 
regime. Thus, a parton shower description of 
high p T emissions is not reliable. Since NLO 
calculations are indeed valid in this region, one 
wants to recover the NLO predictions for high p-y 
emissions after showering. 

Meaningful events: Matching to parton shower is the 
only way to define events at NLO. Without match¬ 
ing, the weights of the real matrix element and 
the subtraction terms belong to different kinemat¬ 
ics and diverge separately. Due to the matching 
scheme, they are combined and one obtains real 
emission phase space configurations with a finite, 
but not necessarily positive, weight. 

We start our discussion of matching from the quan¬ 
tum density matrix. For a generic 2 —> m process at 


NLO, one can write it in a perturbative expansion 

IP) = JP™^ + + IP® i) +0(a 2 s ) . 

Bom,0(l) Virtual, 0(a s ) Real.Otffj) 

Note that we normalize the leading order contribution in 
the counting of the coupling to be of order 1. \f>m) and 
If,®i) correspond to tree level matrix elements, whereas 
|pm ) ) to the one-loop amplitude. The definitions of these 
densities are analogous to the definition given in Eq. (1). 
Based on this quantum density matrix, the observable F 
after showering naively reads 

cr[F] PS =(F\U(t F ,t 0 )\p) = 

°° 

2^7 I [^](^)(®tlt/fe,?o)lp), 

A—m 

where we use the shorthand <!> | = {p, f, s', c', s, c},]. The 
quantum density | p) accounts for the hard matrix ele¬ 
ments for A = m,m + 1. Finally U(tp,to) describes 
the parton evolution during showering, as defined in Eq. 
( 6 ). 

This naive description of the cross section suffers 
from double counting, which is demonstrated as fol¬ 
lows. The evolution equation of the parton shower has 
an iterative solution. Since, at first, we are only inter¬ 
ested in corrections up to 0(a s ), it is sufficient to ex¬ 
pand the evolution equation linearly. Evolving the den¬ 
sity state |p) from to to tp yields 

I p(tp)) = U(t F ,to)\p) « 

|p) + f F dr \<H,(j) - Ip®) + 0(a]). 

Jt 0 

As we can see from the unitarity condition 
(1| [ < W/(r) -'V(r)] = 0, the total cross section 

(1 |p(tx)) is conserved. On the other hand, for inclusive 
observables F one does not recover the NLO prediction, 
because in general (F| [7T/(r) - < V(r)] + 0. Even if we 
would not insists on recovering the NLO prediction, 
the result would have to be considered wrong, since it 
would contain real emission contributions twice: from 
the real emission quantum density |p’°| x ), and from its 
parton shower approximation < 7f/(r)|p^). 

This problem is solved by matching. The solution 
is slightly simpler for processes, which can be defined 
without any cuts at the Born level, e.g. pp —» tt or 
pp —> W 1 VV . At the end, we will obtain a color and 
spin correlated version of the original Mc@NLO for¬ 
malism. The additional parton shower contribution can 
be cancelled by including the same term with an oppo¬ 
site sign in the quantum density state |p). Thus, we can 
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avoid double counting by introducing a modified quan¬ 
tum density state 


Ip) = Ip) - r dT YH,(t) - 'V(r)] \p%>) + 0(a 2 s ) . (7) 
Jt 0 

First notice that (l|p) = (l|p) = cr NL0 is unchanged, due 
to the unitarity condition. On the other hand, consid¬ 
ering U(tf, fo)lp) and expanding the evolution equation 
again shows that the undesired parton shower contribu¬ 
tions are cancelled up to 0(a s ). Notice that this cancel¬ 
lation is non trivial, since the modified quantum density 
|p), depends now explicitly on the parton shower split¬ 
ting kernels and the choice of to. 

Let us investigate the expectation value for an in¬ 
frared safe observable F given by the density state in 
Eq. (7) 




o-[F] =— [dd> m ](F\U(t F ,t 0 )\® m ) 

ml 


(fcmlp!?) + + f F dT(® m \V(T)\p\ 

Jto 

J 


(ill + 1) 

x 


to 

[dO 

m +1 ](F\U(t F , fo)|d> 

m+ 1 ) 
•tp 


IP^r) - f F ^(O m+1 \Hj(t)\ P ^) 

Jto 


Written in this way, the matched cross section suffers 
from infrared divergences in the virtual |p^) and real 
contributions IpJ^jj), which appear in two separate inte¬ 
grals. In the Mc@NLO approach one uses the parton 
shower splitting kernels as subtraction terms, thus one 
drops the infrared cutoff, which is imposed in the parton 
shower, and takes the limit t F —» oo. In the subtracted 
real cross section we can make use of the definition of 
the real splitting operator in Eq. (5) and write 


dr 'Hi(t) — 

/-*oo 

Vs/ I dr 6 (t - ti)®(T - to) = 

\ J o 

y, S/0(p - to). 

I 

Here the sum runs over all external legs and S/ is the to¬ 
tal splitting kernel for a given external leg. We want to 
emphasize that S/ also contains completely finite con¬ 
tributions like the massive g —> QQ splitting, t/ is the 
shower time defined in Section 4.4. Hence 0(h - to) 
represents the ordering of the emissions. The to depen¬ 
dence provides a dynamical restriction of the subtrac¬ 
tion phase space. The real subtracted cross section is 



now finite in d — 4 dimensions, since t/ is allowed to ap¬ 
proach infinity and therefore the subtractions terms can 
resemble the singular limits of the QCD matrix element. 

Integrating the virtual operator 'V(r) without an in¬ 
frared cutoff is more complex, since there is an explicit 
integration over the splitting variables. Thus, we have 
to integrate this part in d = 4 - 2e dimensions to extract 
the 1/e 2 and 1/e poles analytically. 'V(r) takes the form 



dtV(r) = 



dYi S/0(f/ - to) - 


I(to) + K(t 0 ), 


where T/ is the phase space of the additional parton. The 
decomposition of the integrated 'V(r) into I(fo), which 
contains all integrated final-state splittings, and K(fo) 
which contains the initial-state splittings is arbitrary. 
We emphasize this structure only to show that the par- 
ton shower naturally incorporates a subtraction scheme 
as in the Catani-Seymour framework [43]. The only 
part which cannot be generated by the parton shower are 
the collinear counterterms, denoted by P, needed for the 
PDF renormalization. Thus, the matched cross section 
reads 

<r[F] = ( [ ^^(F\U(t F ,to)W m ) 

J ml 

x [(O m |p®) + (Orf) + (® m |[I(f 0 ) + K(f 0 ) + P]|p,?] 


+ i [ ^^(F\u(t F ,t 0 m n+l ) 


So n + 1 )! 

(4WI/0 - |S/|p®)0(f/ - t 0 ) 


( 8 ) 


In practice Eq. (8) is not solved in a single step. Let us 
define the shorthands (as in [4]) 

(3>,„|S) = (® ffl |pL 0) ) + (fcmlpS?) 

+ (ch ffl |[I(fo) + K(fo) + P]|pL 0) ), 


(<S>m + i\H) = (Owlp^) - 2(^ + ilS/lpL 0) )©(f/ - to) ■ 

l 

One can use the total cross section 

-vl 0[1] = if [t /O m ](l|0,„)(O m |5) 
ml J 

+ , * f[dd> m+1 ](W m+l )(® m+ ilH ), 

On + 1)! J 

to generate the hard events according to (<X>„,|S) and 
(0,„ + i|//). The generated events are subsequently in¬ 
terfaced to the parton shower and after the evolution 
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one can investigate the desired observables. Thus, af¬ 
ter showering one has performed the following integrals 

Cr[F] PS = ^7 f [rf® ffl ](F|£/(f F ,/b)l®m)(®m|S) 
ml J 

+ ;-77T f [dO m+ i](FIU(tF,to)l^m+l)(^m+lin) ■ 

(m + 1)! J 

In the case of processes, which require cuts already 
at Born level to yield a finite cross section, one has to 
modify the matching prescription Eq. (8). Naively, one 
would simply introduce a generation cut function given 
by a state (F r \ as follows 

(®»|S) -> (OJSX^Icb,,,) , 

(GW I H) -» (G> m +1 \H)(F /|G> 

m+ 1) • 

Applying the parton shower to these ensembles, shows 
that double counting is still present [72], It turns out 
that it is also necessary to modify (G> m+ i|//) to be de¬ 
fined with account of the generation cut in the subtrac¬ 
tion phase space 


m, — 173.5 GeV, while the bottom quark is considered 
as massless. We use the NLO MSTW2008 PDF set [70] 
with five active flavors and the corresponding two-loop 
running of the strong coupling. We set the renormal¬ 
ization and factorization scale to the top quark mass, 
Pr = Pf = P = m t , and the starting shower time to to 
Tq, with 


e T ° = min 
i*j 


2Pi ■ Pj 

Ql 


where Q ? is the partonic center-of-mass energy. Since 
the ttj process is divergent already at leading order, we 
have to impose cuts on the hard jet in the event gen¬ 
eration. These cuts have to be as minimal as possible 
to ensure the inclusiveness of the events before they 
are passed to the parton shower. We require the recon¬ 
structed jets to have 


Pr(j) > 10 GeV , \y(j)\ < 5 , 

in the event generation, and 

Pr(j) > 50 GeV , \y(j)\ < 5 , 


(G> m+ i| H) -» (Gwlfl) s (G> m+1 |p^ 1 ) 

- ^(G. m+1 |S / |p, ( ° ) )(E / |G / |G), )I+1 )0(r / - t Q ), 

i 

where we have introduced the momentum mapping op¬ 
erator Qi with 

G/|G> m +i) = |G* m ({ P , f}m+l)) • 

as the inverse transformation to that given by the oper¬ 
ator Ri in Eq. (3). It can be shown that this matching 
prescription is correct as long as the generation cuts are 
looser than those implied by the final observable, and 
the latter also only amounts to cuts. 

4.7. Application: ttj production at the LHC 

The matching scheme of the previous subsection has 
been implemented in the framework of Helac-NLO. In 
order to test the implementation, we have chosen to 
study the process of top-quark pair production in asso¬ 
ciation with an additional jet at the Large Hadron Col¬ 
lider. The results reported in this subsection are taken 
from [72], We point out that the NLO QCD correc¬ 
tions to ttj production have been previously obtained 
in Refs. [80-83]. Furthermore, NLO + parton shower 
predictions have been studied in Ref. [84, 85]. 

Results for ttj production are given for pp collisions 
at the LHC with a center-of-mass energy of 8 TeV. The 
top quark is assumed to be stable and its mass is set to 


in the final analysis of several observables. Jets are clus¬ 
tered using the anti -kj jet algorithm [50], with R — 1 
used at both the generation and analysis levels. Only 
particles with pseudo-rapidity I 77 I < 5 are passed to the 
jet algorithm. In the parton shower, the top quark is kept 
as a stable particle (i.e. no decay allowed), hadroniza- 
tion and multiple interactions are not included. In or¬ 
der to address the theory uncertainties we investigated 
the scale dependence on the unphysical scales p and 
Tq. Here p is varied between p = m,/2 and p = 2m,. 
Whereas the parton shower starting time to is varied be¬ 
tween to = Tq/ a/ 2 and to = Tq y/2. 

First we check the total cross section for ttj produc¬ 
tion. Including scale variation of pn we obtain our final 
prediction as 

cr NL0 (pp -> ttj + X) = 86.04+fj 1 ”/^^ pb , 

cr NL0+FS (pp -> ttj + X) = 85 -94^j j 8 y^j %) pb . 

These results are fully consistent, which proves that the 
initial cuts during the generation phase were chosen ap¬ 
propriately. Additionally we observe that the parton 
shower does not improve the theory uncertainty of the 
total cross section at all. This is expected because we 
removed all parton shower effects for the total cross sec¬ 
tion by construction. Nevertheless, we expect some re¬ 
duction of the scale dependence for differential distri¬ 
butions, because of the summation of the leading loga¬ 
rithms. 
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Several distributions showing a comparison be¬ 
tween the fixed order NLO calculation and several 
NLO+PS predictions (aMc@NLO+Pythia8 [20] and 
Powheg+Pythia8 [84]) are shown in Figures 5 and 6. 
pr(j\ - rel) in Figure 6 is the scalar sum of the relative 
transverse momenta of the particles in the first jet, de¬ 
fined with respect to the jet axis in the frame where the 
first jet has zero rapidity 


Pr(j\ ~ rel) = ^ 

iey'i 


I kj x pQj)| 
\p(ji)\ 


where kj is the momentum of the i th particle in /j. 

We observe that distributions shown in Figure 5, 
which are not expected to be affected by showering ef¬ 
fects have, to a good approximation, the same shape as 
at fixed order. The resummation effect is visible first and 
foremost in the transverse momentum distribution of the 
ttj system and other distributions shown in Figure 6. A 
more detailed analysis is presented in [72]. 


5. Conclusions 

We have discussed recent progress towards matching 
next-to-leading QCD calculations for LHC processes 
with parton showers at next-to-leading logarithmic ac¬ 
curacy. Matched NLO+NLL calculations will provide 
accurate predictions for differential distributions and ex¬ 
clusive observables with experimental cuts, and are thus 
essential to fully exploit the potential of the upcoming 
LHC run. 

To facilitate the matching of NLO calculations with 
parton showers, the subtraction terms needed to com¬ 
bine virtual and real corrections should by constructed 
from the splitting functions that define the parton 
shower. We have presented a subtraction scheme based 
on a parton shower with quantum interference and its 
implementation into the Helac-NLO software. The new 
subtraction scheme has been applied to a number of 
challenging processes, including the production of up 
to four heavy quarks at the LHC. The new scheme per¬ 
forms well compared to established methods. It not only 
provides an important internal check of multi-parton 
NLO calculations, but also forms the basis of current 
work on parton shower matching at NLL accuracy. 

First results of such a new NLO+NLL calculation 
have been presented. We find that the resummation is 
important for a wide range of phenomenologically rel¬ 
evant distributions. Work is in progress to further im¬ 
prove on the accuracy of the calculation by adding sub¬ 
leading color effects and spin correlations. 





Figure 5: Differential distributions for ttj production. Comparison 
between NLO and several NLO+PS predictions. 
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